tic;
%------------------------------M, r, y
[nrowAeq, ncolAeq] = size(nAeq);

% check Aeq has full row rank
if rank(nAeq)==nrowAeq
    fullRowRank=true;
end

if fullRowRank
    
%     h = waitbar(0, 'Please wait, geting the solution now...');
    
    % generate a random n x n non-singular matrix M
    M = randi([0,5],ncolAeq);
    
    while det(M)==0
        M = randi([0,5],ncolAeq);
    end
    
    % generate a random n x 1 vector r
    r = randi([0,5],ncolAeq,1);
    
    while nbeq + nAeq * r == 0
        r = randi([0,5],ncolAeq,1);
    end
    
    % generate a random integer between 1 and 50
    gamma = round(rand(1)*49)+1;
    
    cP = gamma * M' * nC;           % c' = gM'c
    
    % -----------------------Q, A', b'
    % generate a m x m non-singular matrix Q (deteminant =/= 0)
    Q = randi([0,5],nrowAeq);
    
    while det(Q)==0
        Q = randi([0,5],nrowAeq);
    end
    
    AeqP = Q * nAeq * M;                    % A' = QAM
    beqP = Q * (nbeq + (nAeq * r));         % b' = Q(b + Ar)
        
    %-----------------------------B'
    [nrowAeqP, ncolAeqP] = size(AeqP);
    nonZero = nrowAeqP;         % find the non-zero of nbeq from end of the vector
    allZero = true;             % determine nbeq is all zero
    
    % find the non-zero of nbeq from end of the vector
    while nonZero ~= 0
        if beqP(nonZero) ~= 0
            allZero = false;            % found & exit loop
            break;
        end
        nonZero = nonZero - 1;
    end
    
    % generate a random n x m lambda
    Lambda =randi([0,5],ncolAeqP,nrowAeqP);
    sum = 0;      % sum of each row
    
    Br = nA * r;              % Br
    % change content of Lambda in order to cause Lambda * b' = Br
    if allZero == false
        if nonZero == 1
            % change first column of Lambda to zero
            for ii=1:1:ncolAeqP                  % row
                if nonZero == 1
                    Lambda(ii, 1) = 0;
                end
            end
        else
            for ii=1:1:ncolAeqP                  % row
                for jj=1:1:(nonZero-1)           % col
                    sum = sum + (Lambda(ii,jj) * beqP(jj));
                end
                Lambda(ii, nonZero) = (Br(ii) - sum) / beqP(nonZero);
                sum = 0;
            end
        end
    end    
 
    AP = (nA * M) - (Lambda * AeqP);         % B' = BM - (Lambda)QAM => BM-LA' =>(A' = AeqP)
    bP = (nA * r) - (Lambda * beqP);         % 0 = Br - Lambda(Q(b+Ar)) =>(Q(b+Ar) = beqP)

    % display transform Time
    toc
    
%     options = mskoptimset('MaxIter',4000000, 'Simplex','on','Write', 'prob.opf','Diagnostics','on');
    options = mskoptimset('MaxIter',4000000, 'Write', 'prob.opf');
    
    tic;
%     [y,fval,exitflag] = linprog(cP,AP,nb,AeqP,beqP,[],[],[],options);
    [y,fval,exitflag] = linprog(cP,AP,bP,AeqP,beqP,[],[],[],options);

    toc
        
    % x2 = round (((M * y) - r) * 100 ) / 100;
    x = (M * y) - r;
    
    n=1;
    nX=[];
    temp=0;
    
    for ii=1:1:cRow
        if n <= length(cVariable)
            if ii == cVariable(n)
%                 if x(ii) < 0.0001
%                     temp = 0;
%                 else
                temp = (x(ii)*100)/100;
%                 end
                nX{n,1} = ['x' num2str(ii) ' = ' num2str(temp)];
                n = n + 1;
            end
        end
    end
    nX
end